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ABSTRACT 

In this paper, we develop a new method for magnetohydrodynamics (MHD) using 
smoothed particle hydrodynamics (SPH). To describe MHD shocks accurately, the 
Godunov method is applied to SPH instead of artificial dissipation terms. In the in- 
teraction between particles, we solve a non-linear Riemann problem with magnetic 
pressure for compressive waves and apply the method of characteristics for Alfven 
waves. An extensive series of MHD test calculations is performed. In all test calcula- 
tions, we compare the results of our SPH code with those of a finite-volume method 
with an approximate Riemann solver, and confirm excellent agreement. 
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1 INTRODUCTION 

It is well known that magnetic fields play an important role in various astrophysical phenomena, such as the formation of stars 
and planets, and high energy astrophysics. Gas dynamics with a magnetic field is well described by magnetohydrodynamics 
(MHD). The MHD equations are too complicated to solve analytically except for special cases, therefore to understand physical 
phenomena involving magnetic fields, numerical simulations are indispensable and powerful tools. A numerical technique to 

solve the MHD equations has been developed successfully in the form of the finite- v olume methods by many authors. 

Smoothed particle hydrodynamics (SPH) is a fully Lagrangian particle method (|Lucvlll977l ; lGingold fc Monaghanlll977l ). 
This Lagrangian nature has major advantages in problems that have a large dynamic range in spatial scale, such as the 
formation of large scale structures, galaxies, stars, and planets. Several authors have tried to apply the S PH method to MHD 
problem s. In this paper, we call SPH for MHD "smoothed particle magnetohydrodynamics" (SPMHD). iPrice fc Monaghan 
( 2004al lbh ha ve developed an on e-dimensional SPMHD scheme. To describe shock waves, they used artificial dissipation terms 
proposed by iMonaghanl (|l997l ) based on an analogy with Riemann solutions of compressible gas dynamics. As the signal 
velocity, they adopted the speed of the fast wave, which is analogous to the sound wave in hydrodynamics (HD). Their 
SPMHD has been sho wn to give good results o n a wide range of standard one-dimensional problems used in recent finite- 
volume MHD schemes. Price fc Monaghan! l|2005l . hereafter PM05) ha ve de veloped a multi-dim ensional SPMHD scheme based 
on the above one-dimensional one. Alternatively. If30rve et alj (|200ll ) and lB0rve et al. (2006) have imple mented an SPMHD 
schem e using a regularization of the underlying particle distribution an d artificial viscosity. Recently. iDolag fc Stasvszvn 
( 20091 ) have implemented MHD in the cosmological SPH code GADGET JSpringel et~aj]|200ll: rSpringellbood) m ainly based 
on PM 05. Broad discussions of the SPH and SPMHD are found in reviews bv lMonagha rJ l|l992l ) JS pringell (|2010al ). and lPrice! 
1 20ld ). 



The SPMHD in PM05 can capture fast shocks accurately. However, since they use the fast wave speed as the signal 
velocity in the artificial viscosity and resistivity terms, Alfven waves become dissipative. In the finite-volume method, it is 
well known that Alfven waves cannot be descri bed accurately if the cha racteristics of Alfven waves are not taken into account 
in the caluculation of the numerical flux (e.g., Stone fc Normar] 19921 ). This is also the case in SPMHD. In this paper, we 



apply the Godunov m ethod to SPMHD. The Godunov method was originally developed in the finite- volume method (jGodunov 
1959; van Leerlll979l ). Unlike artificial viscosity, the Godunov method can, in principle, take into account the minimum and 



sufficient amount of dissipation without any free parameters. In HD, the application of the Godunov method to SPH has 
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been carrried out by Inutsuka ( 20021 . hereafter 102). In MHD, we can also consider the general Riemann problem (RP) with 
arbitrary directions of velocity and magnetic field in both sides. However, it is computationally ex pensive and complex to solve 
the RP because the MHD equations have seven characteristics and non-hy perbolicity. Recently, iGaburoy fc Nitadoril (|201ll ) 



Li 



200^ . Besides the SIM! 



have ap proximate Riemann solver HLLC into a variant particle method i Toro fc Spruce! 19941 ; 

method, Pakmor et al. ( 2011 ) implemented MH D with the appro ximate Riemann solver HLLD ( Mivoshi fc Kusanol 20051 ') in 
unstructured, movin g-mesh code AREPO cod e (|SpringeJl20r0bh . In this paper, we use a simplified approach proposed for 
the finite- volume method by Sano et al. ( 19991 ). In this method, compressible and incompressible parts of the MHD equations 
are completely divided. The former is calculated by a non-line ar Riemann solver with magnetic pressure, and the latter is 
calculated by the method of characteristics (MOC) proposed bv lStone fc Normanl (|1992T ). 

The paper is organized in the following ways. In Section [2] we derive the SPMHD equations from the basic equations 
of MHD. In Section [3] we describe the implementation of the derived SPMHD equations. Various test calculations are 
demonstrated in Section U This paper is summarized in Section [S] 



2 SPMHD EQUATIONS 

The basic equations of MHD can be written as 



d_ 

Tit 





f «" \ 








rp/J.V 







(1) 



where T^" is the stress tensor, 
the specific total energy is given by 

e — P/[(j — l)p] is the specific internal energy, d/dt = d/dt + v ■ V is the Lagrangian time derivative, and we have chosen 
units so that factor of hq does not appear in the equations, where /io is the magnetic permeability. 
In the SPH method, density field is expressed as 

p(x) — rrijW(x — Xj, h(x)), (5) 

3 

where subscripts denote particle labels, rrij is the mass of the j-th particle, W (as, h) is a kernel function, and h is the smoothing 
length that is assumed to depend on x. There are many choices of kernel functions. In the Godunov SPH (GSPH) schemes, 
102 adopted the following Gaussian kernel, 

W M =(^)'e-"* (6) 

where d is the number of dimension. This is because the Gaussian kernel makes the formulation of the GSPH simpler than 
other kernel functions. In addition, the results with Gaussian kernel seem to be better than those with the cubic spline kernel 
in test calculations shown in section [4] In this paper, we follow the choice in 102. 

2.1 Equation of motion 

We define the equation for the evolution of the particle positions as 

Xi = [d 3 x^Wi(x), (7) 



where Wi(x) = W[x — Xi, h(x)]. Substituting equation {TJ into equation 0, one obtains 

= / d 3 x- (V"T F ) Wt(x) = / /xV ( T ^ElM\ _ f d\T^V" (YMp£\ , (8) 
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where we integrate by part. The first term on the right-hand side in equation JSJ becomes the surface integral by the Gauss 
theorem, and vanishes when \x\ — > oo. With equation J5), the last factor of equation |[8j becomes 



) = 7^ mj ' (^'^W^W - WiCaOVWXaO) . (9) 



(10) 



y v f Wi(x) 

Therefore, we can obtain the equation for the evolution of the particle positions, 

•r - - X> ™«.) - wwv^w) ■ 

2.2 Total energy equation 

We define the equation for the evolution of the particle energy as 

E t = J d 3 x^Wi(x). (11) 
Substituting equation |T|) into equation (|11|) . one obtains 



= y ef a:^ (V M T"%") Wi(a:) = - J dTxT^V [^y 1 j ■ (12) 
Using equation ([9]), we can obtain the energy equation of SPH particles, 

d 3 x — {VWi(x)Wj(x) - Wi{x)V"Wj{x)} (13) 

j 

2.3 Induction equation 

The i-th SPH particle is assigned its own magnetic field, Bi. The equation for the evolution of the magnetic field is defined as 
d I B\ f 3 d ( B 



*\ p j-j dx dt\7) Wi{x) - (14) 

Substituting equation (|2| into equation !)14|) . one obtains 

*{t) r I * m ^ r '* Wi{9) (15) 

The right-hand side of equation (|15p can be transformed into the following expression: 

d 3 x— V"B"Wi(x) = / cfx - 1 V ^a)- tPxv»-?-=-Wi(x) (16) 
P J P J P 



We approximate the last term on the right-hand side of equation (|15[1 as follows 



t/ 1 — — Wi(as)d J x = if / — — W i (x)d i x + 0{hT). (17) 

Using equations (|15|l and (117|l and integrating by parts, one obtains the following equation for the evolution of the magnetic 
field: 

s(t) =-S ro i/^^(^-*n{V^(»)Wi(»)-Wi(x)V w Wi(x)}. (18) 



3 IMPLEMENTATION 
3.1 Convolution 

We define s-axis as being along the vector Xi — Xj. The unit vector in the s-direction is TL — (33?' Xj ) j \x% Xj I. The distance 
between i- and j-th particles is Asij = s, — Sj where Si and Sj are the coordinates of the i- and j'-th particles on the s-axis, 
respectively. We need to perform the volume integral in equations (|10|l . (|13[) . and (|18|l . If the smoothing length is spatially 



constant, the volume integral can be done analytically by interpolating physical variables (|Inutsuka 2002). In a similar way, 
one can derive the GSPMHD equations as follows: 
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£ 

dt 
where 



OSi 



(B")v{K)*-in 



(19) 



(20) 



(T MI/ )*, v* , and B* are the values at s = sjy (see Inutsuka 20021 ). and x* = + aSjAt/2 is the time-centered velocity of the 
i-th particle. The quantity Vy(/l) is obtained by the following integration 



p'VifsJ^ta!)/! = Vij(h)W(. 



,V2h), 



(21) 



where 102 used linear and cubic interpolations of p^ 1 (x). The detailed expression of V 2 j(h) is found in 102. However, for 
the case with variable smoothing length, the volume integral cannot be easily performed analytically. 102 proposed a simple 
approximation in which he used hi for the half of the integration space that includes Xi and hj for the other half. In this case, 



dW(A Slj ,V2hi 



+ V%{h j ) 



dW(Asij,V2h:j 



(22) 



This formulation can capture contact discontinuities accurately 1 Cha et al. 2010l : Murante et al. 2011). In this paper, just for 
simplicity, we use the following crude approximation in the volume integral of equations (|10[) . (|13|) . and (|18[) : 



VWi(x)Wj(x) - f.Vffjfi) ~ VW l (x)S(x - Xj) - 5{x - x t )VWj(x) 
Using this, one can get 



dW{As ih hij) 
dsi 



(23) 



(24) 



where fty is an average of hi and hj. This paper adopts fty = (hi + hj)/2. 



3.2 The usage of the Riemann solver 

The equation (|19|l do not include a dissipative process, which is required to describe shock waves. The Godunov method uses 
the exact Riemann solver to include the minimum and sufficient amount of dissipation into the scheme. In the finite- volume 
method, the result of the Riemann problem at cell interfaces is used in the calculation of numerical flux. In the GSPH in 102, 
the values P* and (Pi))* in equations (58) and (59) in his paper are replaced by the results of the RP between the i-th and 
the j-th particles. In the same way, (T^ v )* and (T^ v v")* are replaced by the results of the RP between the i-th and the j-th 
particles. In equation (|19[) . the projection of T M " on the s-axis is found as follows: 



P t -?±\n» + B^, (25) 

where P t = P + Bj_/2, and component parallel (perpendicular) to n is represented by using the subscript of || (_L). The first 
term on the right-hand side of equation (|25|) represents the compressive term working alone the s-axis. In contrast, the second 
term represents the incompressible term working in the perpendicular direction. In the compressible part, we use the result 
of the non-linear RP without Bu , which contains fast shocks, fast rarefaction waves, and one contact discontinuity. From the 
RP, one can obtain P* and The detailed description is shown in Appendix [XI In the incompressible term, MOC is used 
I Stone fc Normanlll992r i. jrom the MOC, one can obtain B* ± and v* ± . The detailed description is shown in Appendix [Bl 



In the calculation of the RP and the MOC, the initial values on each side at s*j are required. In this paper, we adopt 
s*j — for simplicity. It is confirmed that the value of s*j does not affect the results. To make a spatially second-order method, 
we consider the piecewise linear distribution of the physical variables. Using the gradients, the initial values on each side of a 
one-dimensional RP are the average values of each domain of dependence: 

Un = Ut - \ (^) [A«« - dAt] (26) 

Uh = Uj + l(^) j [ASlJ ~ C]At] (27) 

where U = (p, P, v, B), and C is a characteristic speed. In the compressible RP, the speed of the fast wave C = 
y (7P + B 2 )/ p is adopted. In the incompressible RP, the speed of the Alfven wave C = \Bu\/^Jp is adopted. 
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In finite-volume methods with higher spatial accuracy, we need to impose a monotonicity constraint on the gradients of 
the physical variable to obtain a stable description of discontinuities. This is the case in the GSPMHD scheme. The detailed 
description of the monotonicity constraint adopted in this paper is shown in Appendix [Cl 

In actual calculations, we solve the following equation: 

/ -(( J Pt)Rp-P B ,[|)n + B|f(Bx)MOC \ 

~t \ ^ | _ Z_^ m i F n -((-Pt)RP-P B ,||) («||)RP + By(B±)MOC- («i)moc , (28) 

\ B \\ {( u ||)r.p™+ ("i)moc / 

where subscripts "RP" and "MOC" indicate values evaluated using the RP and MOC, respectively, and Pb,|| = {&\ % + B \ j) /4 
and B\ = (B ht +B hj ) /2. 



d 




3.3 Variable smoothing length 

In this paper, the variable smoothing length is used to obtain large dynamic ranges. The smoothing length of the i-th particle 
is determined iteratively by 

/ \ 1/d 

h * = Ch WT rr , (29) 

where Ch is a parameter. The density of the i-th particle is evaluated by 

pi = uijW^Xi - x 3 , hi). (30) 

3 

The Gaussian kernel is not truncated at a finite radius but has infinite range. In practical calculations, we ignore the contri- 
bution from the j'-th to i-th particles if \xi — Xj\ > 3.1h because exp(— 3.1 2 ) ~ 6.7 x 10~ 5 is sufficiently small. The number 
of neighbours becomes ~ 6Ch in ID, ~ 30Cf in 2D, and ~ 124C^ in 3D schemes. In this paper, we present the results of 2D 
test calculations and adopt Ch = 1.2, indicating that the average neighbour number is ~ 43. 



3.4 Corrections for Avoiding Tensile Instability due to Magnetic Force 



From equation (|25[) . one can see that the stress tensor can be negative when the plasma beta ft = 2P/B 2 is low . This causes 
unph ysical clumping of SPH particles ( Monaehanll 19921 ) . This numerical instability is called "tensile instability" I Swede et al. 
1995). The tensile instability arises from the fact that the SPH expression of V • B is not completely zero. In finite-volume 



methods, it is well known that nonzero V • B pro duces an unphysical force along the magnetic field and causes large errors 
in the simulations when using conserv ative form jBrackbill fc Barned[l980h . Many authors proposed methods for v anishing 
V • B, such as the projection method ( Brackbill fc Barnes 1980l ). the constrained transport I Evans fc Hawlev 19881 ). and so 



In SPMHD, V • B inevitably has some amount of numeric al noise that causes the tensile instability to occur even if 
divergence-cleaning methods are used in the con servation form JPrice fc Monagharj[2005l ). Therefore, several methods have 
been proposed to suppress the tensile instability. Phillips fc Monaghan ~(l7985h proposed t hat the stress tensor make positive 
by sub tracting an constant value from the stress tensor (also see Price fc Monaghan 20051 ) . As another approach, Borve et al 



I 200lh suggested that the monopole source terms are explicitly subtracted from the equation of motion and the energy equation 
as follows: 
dv» 1. 
~dt 

and 



P 



V B, 



(31) 



dt p 



P 



(32) 



The left-hand side of equation (|31[) corresponds to the Lorentz force ((V x B) x B)/p. This source terms signific a ntly s tabilize 



the tensile instability. This formulation is the same as so-called 8- wave formulations proposed bv lPowell et all (|1999T ) in the 



finite- volume method. This formulation is numerically stable, and the value of V-i? keeps zero within truncation error without 
any cleaning methods for simple test problems. However, in realistic 3D problems, satisfaction of the divergence constraint is 
not guaranteed. Thus, we need to monitor the value of V ■ B in actual calculations. To derive the GSPMHD equations, we 
take the convolution of the source term in equation (|31[) . 



■BV ■ BWi(x)d 3 x ~ -Bi 



Bi / B V 



Wi(x) 



d 3 x. 



(33) 
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Using equation ©, one can get 

= £>i {( T ^yn - B^B;}F ij: (34) 

3 

In a similar way, the energy equation is also obtains as 

E$ = J2™i {(T^v")* - Bi ■ xir } / .,. (35) 

3 

The major disadvantage of this correction is violation of the momentum and the total energy conservation. Actually, in all 
SPMHD schemes without the tensile instability, the conservations of the momentum and/or total energy are sacrificed. 

3.5 Divergence Error Estimate 

In SPMHD, there are many choices of V • B estimate. In this paper, we use the following expression: 

(36) 



(^-V-B^j = J iv ■ BW{x-x l )d A x = ^m 3 BlF, ] . 



3 

In order to estimate errors in the V • B = constraint, we monitor 5b defined as 

*B = ^E ^T B f /P H (37) 

JVtot * — ' \tSi\ 
i 

where JVtot is the tot al number of the SP H particles. This estimate is the same as that in the correction terms in section [3741 
Note that PM05 and Bdrve et al. ( 20061 ) adopted different choice, 



(V ■B) i = — S2m J {B t -B J )-V t W(x l -x J ,h t ). (38) 

pi 

3 

The divergence error estimated in equation (|38[) tends to be smaller than those in equation (|36|) . 



4 NUMERICAL TESTS 

In this section, we show the results of test calculations in the two-dimensional GSPMHD method, starting with the convergence 
test. 

4.1 Convergence Test 

In this section, we check the accuracy of the GSPMHD. A good problem for the convergence test is the propagation test of 
linear MHD waves in a uniform media. First, we define the unperturbed state. A uniform and rest gas is considered (p = 1, 
P = 1, v — 0). The magnetic field is parallel to the ^-direction, and its magnitude is l/y/2. The simulations are performed 
in the square domain, x, y £ [0, \/2] • The particles are uniformly spaced on a cubic lattice with sides parallel to the x- and 
y-axes. A periodic boundary condition is imposed. 

We consider linear MHD waves having a wavenumber of k — 27r(l/v / 2, l/v2, 0), indicating that the simulation domain 
contains two wavelengths. It is well known that MHD linear waves consist of the fast, Alfven, and slow waves with phase 
velocities given by 1.4, 0.5, and 0.46, respectively in this configuration. We add the initial perturbation based on the eigenmode. 
In the fast and slow modes, all fluctuations lie in the {x,y) plane. The amplitude of the density fluctuation is set to 10 -3 . 
In the Alfven wave, only v z and B z fluctuate. The amplitude of the fluctuation of J3 Z is set to 10 -3 . The vector component 
parallel (perpendicular) to the wavenumber is expressed in terms of the subscripts £ (t/>) in the (x, y) plane. 

As a measure of the error, we introduce the error vector defined as 

iVtot 

e= 1 -—y^\U le{ (xi)-Ui(x i )\, (39) 

JVtot ' ' 

i = l 

where U = (p, ti£, v^,, B^, E) for the fast and slow waves, U = (v z , B z ) for the Alfven wave. As a reference solution, J7 re f, 
we adopt the results with JVtot = 512 x 512. To eliminate the error coming from At, At is set to the small value of 3 x 10 -4 in 
all resolutions. The error vector is evaluated after 100 time steps at various resolutions. In Fig. [1] the norm of the error vector 
is plotted as a function of the average smoothing length, which represents the resolution for the fast (the circles), Alfven 
(triangles), and slow waves (boxes). In the scheme having second-order of spatial accuracy, |e| is expected to scale as h 2 . Fig. 
[1] shows that the error is proportional to h 2 for all wave modes. Therefore, it is confirmed that our GSPMHD is spatially a 
second-order scheme. 
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Figure 1. Results of convergence test of the fast (circles), Alfven (triangles), and slow waves (boxes). The abscissa indicates the average 
smoothing length that represents the resolution. The ordinate indicates the norm of the error vector e. The upper and lower solid lines 
represent the lines of <x h~ 2 and of <x h , respectively. 
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Figure 2. Results of the circularly polarized Alfven wave test after 5 periods. The ordinates indicate the magnetic field perpendicular 
to the wave number, k. The abscissas indicate the projected coordinate in the direction of k. The analytic solution is shown by the solid 
line in each panel. All particles are plotted by the circles. The results are shown at three different resolutions TVtot = 16 X 32, 32 X 64, 
and 64 X 128 from bottom to top. The left and right panels correspond to the hexagonal lattice and th random distributions in the initial 
condition, respectively. 



4.2 Non-linear circularly polarized Alfven wave 



Tothl (|200C ) investiga ted a non- linear circularly polarized Alfven wave that is one of the exact solutions of the non- linear MHD. 
Following Toth (2000), we set the following initial condition. The Alfven wave propagates toward an angle a = n/6 with respect 
to the x-axis. The initial condition is p = 1, P — 0.1, B^ = 1, = — 0.1 sin (2-irx^), and v z — B z = 0.1 cos (2nx^), where 
= a; cos a + ysina. In order to investigate the effect of the particle distribution, two kind of initial particle distributions 
are considered. One is the ordered distribution (a hexagonal packed lattice), and the other is the random distribution that is 
relaxed until the density dispersion is sufficiently small. In each particle distribution, we calculate this test at three resolutions, 
16 x 32, 32 x 64, and 64 x 128 particles. 
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Figure 3. Results of the circularly polarized Alfven wave test with (a) the cubic spline kernel and (b) the Gaussian kernel after 1 period. 
The total particle number is 64 X 128. The half of the calculation region (0 ^ ^ 1) is plotted. The solid line incidate the analytic 
solution in each panel. In each panel, the results are shown with different smoothing length Ch = 1.0 (the circles), 1.2 (the triangles), 
and 1.5 (the boxes) (see equation H29I I), 




Figure 4. Results of the circularly polarized Alfven wave test without a monotonicity constraint for iVtot = 64 X 128 (the circles) after 
5 periods. The solid line indicates the analytic solution. 



The results are shown in Fig. Rafter five periods. The abscissa and the ordinate denote the projected coordinate x% and 
the perpendicular magnetic field in the [x, y) plane, respectively. The exact solution is shown by the solid line in each 
panel. All particles are plotted by the circles. The results are shown at three different resolutions AT to t = 16 X 32, 32 x 64, and 
64 x 128 from bottom to top. The left and right panels correspond to the hexagonal lattice and the random distributions as the 
initial condition, respec tively. Fig. [2] also shows t hat the results with different initial particle configurations agree extremely 
well in each resolution. iPrice fc Monag han (2005) performed the same test. In their Fig. 6, the phase error is found even in 
the highest resolution. In all panels of Fig. [5] the phase error is not seen in our GSPMHD. The phase error may come from the 
fact that PM05 used the cubic spline kernel. The results with the cubic spline and the Gaussian kernels after one period are 
shown in Figs. [3^, and[3p for N to t = 64 x 128, respectively. To see the phase error clearly, Figs. [3] are plotted for ^ ^ 1. 
In each panel, we show the results with different smoothing length Ch = 1.0 (the circles), 1.2 (the triangles), and 1.5 (the 
boxes) (see equation (|29[l ). The solid line indcates the analytic solution in each panel. From Fig.[3K, even after one period, one 
can see that the cubic spline kernel gives relatively large phase errors, the values of which depend on the smoothing length, 
Ch- On the other hand, the Gaussian kernel shows sufficiently small phase errors compared with the cubic spline kernel, and 
the phase errors are nearly independent of Ch- Therefore, the Gaussian kernel is superior to the cubic spline kernel in the 
propagation of Alfven waves. 

Fig. [2]shows some errors around the extremal points of at a;^ = 0.25, 0.75, 1.25, and 1.75. This error comes from the 
strange "clipped" shape of the wave compared with the exact solution in Fig. [2] while it was not found in PM05. This is caused 
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by the monotonicity constraint on the gradients of the physical variables, required for a stable description of discontinuities 
(see section [3T2"1) . Because the monotonicity constraint makes the scheme's spatial accuracy first-order around extremal points, 
the profile around extremal points dissipates preferentially and is flattened as seen in Fig. [2] Fig. [4] shows the results without 
the monotonicity constraint for iVtot = 64 x 128 (the circles) after five periods. One can see that no deformation arises. 
Deformation of the wave shape is al so found in fin ite-volume methods using the monotone upstream-centered scheme for 
conservative laws (MUSCL) method ijvan Leerll 19791 ). Apart from the wave shape, the waves in GSPMHD are less dissipated 
than those in PM05, who adopted the artificial resistivity. 



4.3 Shock Tube Problems 



MHD shock tube problems are widely used to test numerical codes. In this section, we calculate two shock tube tests. 

First, we perform a shock tube where initial states are given by (p, P, v x , v y , v z , B y , £> z ) = (1.08, 0.95, 1.2, 0.01, 0.5, 
3.6/V47T, 2/^/47) for x < and (1, 1, 0, 0, 0, 4/V47T, 2/-v/4tt) for x > with B x = 2/V47T. Initially, SPH particles are 
distributed in a hexagonal lattice with the particle separation of 4 x 10 -3 . The particle separation in the x-direction for x > 
widens slightly to obtain the initial density disc ontinuity. The rectangula r do main is [—0.74, 0.51 x [-3.2 x 10~ 2 ,3.2 x 10" 2 ]. 
The same shock tube problem was presented by Dai fc Woodward! (1994) and Rvu fc Jonesl (1995) in finite- volume methods 
and PM05 in SPMHD. The exact solution consists of two fast shocks, two rotational discontinuities, two slow shocks, and one 
contact discontinuity. Fig. [5] shows the results of the GSPMHD at t = 0.2. The solid gray lines indicate the exact solution. 
One can see that GSPMHD describes all discontinuities very well. Our GSPMHD can resolve the rotational discontinuities 
and the slow shocks although they are smeared out in Fig 9 of PM05. This may illustrate the contrast that our GSPMHD 
takes into account the characteristics of Alfven waves and PM05 use an artificial resistivity in the induction equation. 

Next, we perform a shock tube test contains stronger shocks than the previous one. The initial states are given by 
(p,P,v x ,v y ,v z ,By,B z ) = (1,20, 10, 0, 0, 5A /47T, 0) for x < and (1 1 -10, 0, 0, S/y 7 !^, 0) for x > with B x = The 
same shock tube problem was presented bv lDai fe Woodward! l| 19941 ). iRvu fc Jonesl |l995l ) and iTothl (|200oh in finite- volume 
methods. The exact solution consists of two fast shocks, a left-propagating slow rarefaction wave, a right-propagating slow 
shock, and one contact discontinuity. The initial particle distribution is a hexagonal lattice with the average particle separation 
of 5.4 x 10~ 3 in [-1, 1] x [-7.2 x 10" 2 , 7.2 x 10" 2 ]. Fig. [6] shows the results of the GSPMHD at t = 0.06, with the solid 
gray lines indicating the exact solution. This figu re shows tha t GSP MHD can reproduce the exact solution and describes 

The fast shocks can be resolved by small number of 



all discontinuities better than those in PM05 and 



arve et al 



particles. In the finit e-volume method, Toth ( 200ol ) reported a relatively large error of B x in his Fig. 13 where the 



conservative method I Powell et al. 19991 ) is used. However, even in this kind of shock tube problem with strong shocks, our 
scheme shows the error in B x less than 1 percent except for the vicinity of the discontinuities. 



4.4 Orszag-Tang Vortex 



The next test is the Orszag-Tang vortex problem that was originally investigated bv lOrszag fc Tang II197E ) in incompressible 
MHD flows. This problem is a standard two-dimensional test for compressible MHD schemes ( Toth 2000). This calculation 
is performed in [0, 1] x [0, 1] domain. In all boundaries, periodic boundary conditions are imposed. The initial conditions are 
given by p = 25/(36tt), P = 5/(12tt), 



v(x,y) — (sin(27rj/), sin(27ra;), 0) , and B(x,y) — 



(— sin(27ry), sin(47ra;), 0) . 



(40) 



Although the initial velocity and magnetic field are not random, the system moves into turbulence through non-linear inter- 
action of MHD waves. Fig. [7| shows that contour maps of the test at t = 0.5 for (a)density, (b)pressure. (c)ma gnetic energy 
B 2 /2, and (d)specific kinetic energy v 2 /2. Fig. can be directly compared with Fig. 22 of Stone et al. (2008) and one can 
see that the agreement is excellent. To view the results quantitatively, we compare the horizontal cuts of the temperature 
for GSPMHD and a finite-volume method with HLLD Rieman solver and the constraint transport method (provided by Dr. 
T. Matsumoto) in Fig. [8] at (a)y — 0.5, (b)y = 0.427, and (c)y — 0.3125. The black and grey solid lines denotes the results 
with the GSPMHD and the finite- volume method with the resolution of 256x256. One can see that the profiles between the 
two methods show good agreement except for the peak at x = in the y — 0.5 slice. Our scheme does not strictly conserve 
the total energy, E to t = rriiEi. The time evolution of the relative error of E to t is presented in Fig. [9K- One can see that 
the error of E to t is sufficiently small. Fig. [9p shows the time evolution of the divergence error, which is maintained at an 
acceptable level ~ 1 parcent. The distribution of the divergence error localizes at shock fronts. The spatial distribution of the 
SPH particles for /i, |V ■ B\i/\Bi\ > 0.05 is shown in Fig. (9ji. The divergence error is highly localized at discontinuities. This 
error comes from the irregular particle distribution. 
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Figure 5. Results of shock tube problem at t = 0.2. The initial condition is (p, P,v x ,v y , v z , B y , B z ) = (1.08, 0.95, 1.2, 0.01, 0.5, 3.6/V4tt, 
2/V?n) for x < and (1, 1, 0, 0, 0,4/ v^, 2/ v / 4^) for x > with B x = 2/Vtn. The circles indicate results of GSPMHD. The solid gray 
lines indicate the exact solution. 



4.5 Rotor 

The MHD rotor problem was introduced by Bals ara fc Spicei (1999) to test propagation of strong torsional Alfven waves. The 
computation domain is a square unit [—0.5,0.5] x [—0.5,0.5]. This problem consists of a dense and rapidly rotating cylinder 
(rotor) embedded by a rarefied uniform medium. The initial conditions are given by 

p = 10, v = (-2y/r ,2x/r ,0) for r = y/x 2 + y 2 < r , (41) 
p=l + 9/(r), v = f(r)(-2y/r,2x/r,0) for r <r<n, (42) 
and 

p=l, v = (43) 

for r > n, where f(r) — (ri — r)/(n — ro), ro = 0.1, and n = 0.115. The pressure and the magnetic field P = 1, 
-B = (5/v47T, 0, 0) are uniform, and the adiabatic index is 7 = 1.4. All SPH particles are assumed to have the same mass. 
Therefore, the number density of SPH particles in the rotor is larger than that in the ambient gas. The SPH particle mass is 
determined so that the resolution in the ambient gas is as large as 256 x 256, leading that the total particle numb er is 86968. 
The initial particle distribution is constructed by using a relaxation method presented in IWhitworth et all jl995l ). 



Godunov smoothed particle magnetohydrodynamics 11 




Figure 6. Results of shock tube test with the initial condition (p, P, v x , v y , v z , B y , B z ) = (1, 20, 10, 0, 0, 5/v4tt, 0) for x < and (1, 1, 
-10, 0, 0, 5/VItt, 0) for x > with B x = 5/\/4tt at t = 0.06. The circles indicate results of GSPMHD. The solid gray lines indicate the 
exact solution. 



Fig. 1101 shows the contour maps of (a)density, (b)gas pressure, (c) Mach number \v\/ y/ yP/p, and (d)magnetic energy 
B 2 /2 at t = 0.15. The contour levels are the same as those in Fig. 18 o flTothl (120001). From Fig. 10| the results agree with 
Totbl 120001 ) quite well. Compared with other SPMHD schemes, such as brice fc Monaghar] (|200ol ); lB0rve et all 120061 1, the 



contours appear to be smoother. 

To compare with the finite-volume method in more detail, we plot horizontal slices of the rotor problem at y = 0.5 (first 
row) and at x — 0.5 (second row) in Fig. 1111 In each panel, the black and grey lines indicate results with GSPMHD and the 
finite-volume method, respectively. One can see that the agreement between the two methods is quantitatively excellent in all 
variables. Since the GSPMHD is a Lagrangian method, the resolution is better in the dense ring while the GSPMHD gives 
more diffusive results in the narrow region with low density ahead of the dense ring . Fig. [12] is the same as Fig. [9] but for the 
rotor test. In this test, the energy and the divergence error is maintained in a sufficiently low level. 



4.6 Blast Wave in a Strongly Magnetized Gas 

The next test is blast waves that propagates into a strongly magnetized gas (jBalsara 



Spicer 19991 ; Londrillo fc Del Zannal 
2000). The calculation region is a square domain of [—0.5, 0.5] x [—0.5, 0.5]. In this problem, an overpressured hot region with 
pressure of Phot is set within r < ro, where r — \J x 2 + y 2 . Around the hot central region there is a rarefied ambient gas 
with a pressure of P a mb- The density p = 1 is spatially uniform. The initial uniform magnetic field B = Bq(1/\/2, l/y/2, 0) 
makes the angle n/4 with the z-axis. Here, ro, Phot, Pamb, and Bo are parameters. We consider two cases: moderate- /3 case 
and \ow-/3 case. 



4-6.1 Moderate-ft Case 



We adopts ro = 0.1, Phot = 10, Pamb = 0-1, and Bp = 1. Therefore, the plasma /3 of the ambient gas is as low as 0.2. These 
parameters were adopted in iGardiner fc Stone) |2005h . Fig. 1131 shows that the contour maps of the blast wave at t = 0.15 
for (a)p, (b)P, (c)i> 2 /2, and (d)B 2 /2. The total particle number is as large as 256 x 256. One can see the shock structures 
around the elongated hot bubble along Bo- The fast (slow) shock propagates tow ard the direction parallel (perpendicular) 
to Bq- This figure can be directly compared with the bottom column of Fig. 28 in lStone et al.l ( 20081 ) . One can see that the 
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Figure 7. Contour maps of the Orszag-Tang vortex test at t = 0.5 for (a)density, (b)pressure, (c)magnetic energy B 2 /2, and (d)spccific 
kinetic energy v 2 /2. 

GSPMHD 256x256 




0.0 0.5 1.0 0.0 0.5 1.0 0.0 0.5 1.0 



Figure 8. (a)Horizontal slices of the temperature in Orszag-Tang vortex at t = 0.5 taken at (a)j/ = 0.5, (b)y = 0.427, and (c)y = 0.3125. 
The black and gray lines in each panel denote the results with the GSPMHD and the finite-volume method, respectively. 
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Figure 9. Time evolution of the error of (a)the total energy | Stat /-Etot (t = 0) — 1| and (b)the divergence error 5b for the Orszag-Tang 
Vortex test. (c)The spatial distribution of the SPH particles for /i;|V ■ > 0.05 




Figure 10. Contour maps of rotor problem at t = 0.15 for (a)density, (b)gas pressure, (c)Mach number, \v\/ \J 7P '/ p, and (d)magnetic 
energy B 2 /2. The contour lines are the same as those in Toth (2000). 




Figure 13. Contour maps of the moderate /3 case at t = 0.15 for (a)dcnsity, (b)gas pressure, (c)specific kinetic energy (u 2 /2), and 
(d)magnetic energy (B 2 /2). The 30 contour lines arc shown for the ranges 0.14 < p < 2.78, < P < 0.95, < v 2 /2 < 0.37, and 
0.105 < B 2 /2 < 1.4. 
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Figure 14. Horizontal slices of the blast wave in the relatively low f) case at (a)a = 7r/4 (first column) and at (b)a = — tt/4 (second 
column). The density, pressure, velocities, magnetic fields are plotted. In each panel, the black and gray lines indicate the results with 
the GSPMHD and the finite-volume method, respectively. 
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Figure 15. The same as figure[9]but for the blast wave test. 



contour maps is quite similar to those in Stone et al. ( 2008h except for the central rarefied hot bubble where GSPMHD is 
more diffusive owing to its Lagrangian nature. 

For comparison with the finite- volume method in detail, we consider slices of physical variables passing through the centre 
(0, 0). To characterize the direction of the slice, we introduce an angle a that is the angle between the slice and the a;-axis. Fig. 
If 41 shows the results for the cases with a = n/4 and a = — n/4, which correspond to the directions parallel and perpendicular 
to the initial magnetic field, respectively. The subscripts || and _L represent the components parallel and perpendicular to the 
direction of the slice. The black and grey lines indicate the results with GSPMHD and the finite- volume method, respectively. 
One can see that the results with GSPMHD agree with those of the finite-volume method very well except in the central 
region, as mentioned above. In the profile of the parallel magnetic field B\\ for a = 7r/4, there are some wiggles near the 
contact discontinuity. This comes from the pressure jump in the initial condition, and does not serious because B\\ agrees with 
that in the finite- volume method in the other places. If the pressure jump is smoother, the wiggle becomes small. Fig. 1151 is 
the same as Fig. [9] but for the blast wave test. In this test, the energy and the divergence error is maintained in a sufficiently 
low level. 



4.6.2 Low 13 Case 

We adopts ro = 0T25, Phot = 100 , Pamb = T and Bp = 10. Therefore, the plasma 8 of the a mbient gas is as low as 0.02. 
These parameters were adopted in lLondrillo fc Del Zannal ^2000) and iGardiner fc Stond (|2005l ). Fig. 1161 is the same as Fig. 
1131 but for the low B case. One can see a stronger slow shock along the initial magnetic field than that in the previous case. 
Fig. 1171 shows the slices of the physical variables along a = ±7r/4 with respect to the x-axis, and is the same as Fig. 1141 but 
for the low B case. One can see that the results of GSPMHD coincide with those of the finite- volume method very well also 
in the low B case. 



5 SUMMARY 

In this paper, we developed a new SPMHD scheme with the Godunov method. To take into account the physical dissipation, 
we consider the non-linear RP with magnetic pressure and the MOC in the interaction between the SPH particles instead of 
artificial dissipation used in previous works. Using the MUSCL method, the spatial accuracy of our scheme attains second order 
accuracy OQi 2 ) that is confirmed in the convergence test of linear MHD waves (see section |4|. From several test calculations 
in section [4] it is confirmed that our method can capture all MHD discontinuities more accurately than previous proposed 
methods. The GSPMHD can provide results comparable to finite-volume methods with approximate Riemann solvers. We 
will apply the GSPMHD to astrophysical problems where the Lagrangian description has advantages. 

The GSPMHD described in this paper loses strict conservation property with respect to both momentum and energy. 
This is to avoid the tensile instability in strong magnetic fields. Although conservation errors are sufficiently small in the 
test calculations (see section |4|, the non- conservative formulation can be problematic in long-term calculations. Thus, further 
investigations are needed to improve the conservation property together with the better performance in low-/3 plasma cases. 

In the GSPMHD, the Gaussian kernel is used. Other kernel functions (e.g., the cubic spline kernel) can be applied easily 
in our GSPMHD if one use equation (I24|l or other symmetrization of the kernel function. However, in the SPMHD, choice of 
kernel functions may be important in contrast to the HD case. In section [4721 it is shown that the cubic spline kernel brings 
relatively large phase errors into the propagation of Alfven waves. In shock tube tests, we confirm that the results with the 
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Figure 16. The same as Fig. 1131 but for the low /3 case at t = 0.02. The 30 contour lines are shown for the ranges 0.233 < p < 3.31, 
32.1 < P < 1.1, < v 2 /2 < 13, and 24.5 < B 2 /2 < 76. 
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Figure 17. The same as Fig. [Tl]but for the low /3 case at t = 0.02. 



cubic spline kernel are worse compared with the Gaussian kernel. To obtain reasonable results with the cubic spline kernel, 
one need a large neighbours hi > 1.5(mi/pi) 1,/2 in the two-dimensional code as suggested in PM05. Thus, we recommend the 
Gaussian kernel in the GSPMHD. 
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Figure Al. Schematic picture of the Riemann problem. "CD", "FS", and "FR" denote the contact discontinuity, the fast shock, and 
the fast rarefaction wave. 



APPENDIX A: RIEMANN SOLVER WITH MAGNETIC PRESSURE 



In this Appendix, we present the non-linear Riemann solver. Fig. IA1I shows that the schematic picture of the non-linear RP. 
Initially, we consider two uniform states J7l, Ub, that are separated by the discontinuity at m = 0, where U = (p, P, v, B±) 
and m = J" pds is the mass coordinate. The magnetic field in the s-direction, By , is assumed to be zero. The RP depends 
only on B\. Since dv±/dt = for Bu — 0, v± is constant spatially and temporary in each side even if v± has a discontinuity 
at m = 0. Therefore, v± does not affect the RP, suggesting that we can set v± — without loss of generality. 

In this RP, the fast shock (FS) or the fast rarefaction (FR) waves propagate o utwar d as shown in Fig. IAII This config- 
uration is the same as the RP used in the Godunov method of the HD ( Godunov 19591 ). The RP is separated into the left 
and the right intermediate state, U1 and J7 R , by the contact discontinuity (CD). At the CD, since the total pressure and the 
velocity is continuous, one can get the following relations, 



p: = pi + 



= Pr + 



(Bl 



(Al) 
(A2) 



In order to take into account the FR exactly, we need the numerical integration that is computationally expensive. Therefore, 
we treat the FR as "rarefaction shock". This mean that the shock jump condition is used in the case of P* < Pth or PtR. 
This treatment is reasonably accurate because the tangential lines of the Hugoniot curve and the adiabatic curve coincide at 
any point in (p, P t ) plane. 

From equation {TJ, the relations between the jump of the total pressure and velocity across the right-facing fast shock 
and the left-facing fast shock are given by 



P t * - PtR = M R - t>|| R ) , and P* - P tL = -M L (ujj - «|| L ) 



(A3) 



respectively, where Mr and Ml are the Lagrangian speeds of the right-facing and the left-facing shocks, respectively. From 
equation (|A3|) . the total pressure and the velocity in the intermediate state are given by 



Pt* 



1 



1/Mr + 1/Ml 
1 



PtR PtL 

Mr Ml 



IWllR 



[MrD|| R + M L W|| L - (Pi 



11 Mr + Ml 

The Lagrangian shock speeds Ml, Mr are derived from the jump conditions across the shock, 



9 ^ B \ 

Pvl + P + ^t 
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(A6) 



[v\\B x ] = 
Using equations 

pa 

4 



P v \\ lP 
2 7-1 



■vuBi 



and (|A7|I . one can get 

Ml = % [( 7 - 3) P ta + (7 + 3) P; - (7 ~ 2) B 2 ±a + VL~\ 



(A7) 



(A8) 
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Figure Bl. Schematic picture of the method of characteristics for Bu > 



D a = {( 7 + 1) P ta + ( 7 - 1) P:} 2 - (7 ~ 2) Bi a {2 ( 7 - 3) P ta + 2( 7 + 3)P t * - (7 - 2) flL} , (A9) 

where a =L and R. From equation (IA8|) , since Ml and Mr depend on P t *, equation ()A4() is non-linear with respect to P* . 
Therefore, we solve equation (|A4[) iteratively by the following procedure. First, the total pressure P t *^ = (PtR + Pth)/2 is 
inserted into the right hand side of equation (|A4|) . Then, we can get P^' 2 - 1 that is also inserted into equation (1A4[) to get P^' 3 '. 
The iteration is continued until the desired accuracy is reached. Finally, the velocity v?, is obtained from equation (|A5[) . 



APPENDIX B: METHOD OF CHARACTERISTICS 

In this Appendix, MOC is briefly reviewed. We consider the propagation of the Alfven wave along the s-axis. For simplicity, 
we consider the case with B\\ > 0. More general expression including the case with Bu < is presented later. The MHD 
equations for a one-dimensional (the s-direction) incompressible fluid are given by 



dvj_ _ P|| 8Ba 
dt ~ 



and 



dB, 



p ds dt 
equations (|B1[) can be written as 



dv± 
~d7 



dv± 
dt 

and 

dv± 
~dT 



1 dB± 
yfp dt 



+ 



P 



(dv ± 



1 dBj 



VP\ds y/p ds 



= 0, 



+ 



1 dB± 

VP dt 
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1 dB± 

VP 9s 



From equations (|B2|) and (|B3|) . one can see that 



dj+ = dv± — 



dB A 



and dj- — dv± + 



dB, 



(Bl) 



(B2) 



(B3) 



(B4) 



y/P VP 

are constant on a trajectories with ds/dt — Bu/^/p and ds/dt — —B\\/^/~p~, respectively. Fig. IB II shows the schematic picture. 
The partially updated values v*^ and B* ± at t + At/2 can be obtained by extrapolate back in time along C+ and C_ to the 
present time step where all variables are known. The positions s+ and s_ are the foot points of C+ and C_ intersect at s* 
on t + At/2, respectively (see Fig. IB1|) . Using the Riemann invariant dj±, the characteristic equations along C+ and C_ are 
given by 



B", 



= 0, and v± — v , + 



B*, — B , 



= 0, 



respectively. From equations (|B5[) 

'pi 



B* ± and v*j_ are given by 
B- 



v 



■ v 



(B5) 



(B6) 



and 
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Figure CI. Schematic picture of the linear interpolation 



p+vj + Vp~ v 1 ~ B 1+ B a 



(B7) 



respectively. So far, we consider only for the case with Bu > 0. For the case with Bu < 0, the Alfven wave propagates in the 
opposite direction of the s-axis. Therefore, the positions of s + and s~ replace each other. The general expressions are given 
by 



Bl = 
and 
v*± = - 



1 1 
-= + 



+ sgn(B||) (-«!+' 



i+v+ + y/p~=v± + sgn(B||) (-B+ + Bl) 



(B8) 



(B9) 



Vp + + VP' 

where sgn(_B||) is the sign of Bp In actual calculations, as By, we use the following simple average value, (B\\ t i + /2. 



APPENDIX C: MONOTONICITY CONSTRAINT 



To solve the RP and the MOC in the interaction between the i- and j-th particles, we need to evaluate physical variables Ul 
and Ur at s = (sj + Sj) /2 as shown in section \S7Z\ In the GSPMHD, a second order spatial accuracy can be achieved using a 
piecewise linear interpolation of the physical variables to determine Ul and J7r. Fig. I C 1 1 shows the schematic picture of the 
linear interpolation of a physical variable Q that is one of U. The slope of the linear interpolation of Q is simply assigned by 
the gradient at each particle's position that is evaluated as 



(dp) 



m k VW(xi — x k , hi) 



for Q = p, otherwise 



(8Q\ 
\dsj t 



V — (Q k - Qi) VW(xi -x h ,hi 
^ Pk 



K = Qj + -AQj, and Q R = Qi - -AQ, 



Using the gradient, the values at left and right side in the RP are given by 
1 . „ , „ 1 

respectively, where 
AQi = (j^-) As iJ5 and AQ 3 



fdQ\ 



As,. 



(CI) 



(C2) 



(C3) 



(C4) 
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If equation (|C2I) is directly used in the derivation of Ql and Qr, unphysical numerical oscillations arise l|van Leerlll979h 



In order to obtain stabl e second-order scheme, we need to impose a monotonicity constraints on AQ. In the finite- volume 
method, van Leerl (|l979h proposed several monotonicity constraints. We apply one of them to the GSPMHD as follows 



mm{2\Q i -Q j \, \AQi\, 2|AQ;|}sgn(AQ0 
(AQ)-H if sgn(Q l -Q J )=sgn(AQ l )=sgn(AQO, (05) 

otherwise, 

where AQi satisfies 

AQi = ~ Qj) + ^Q'j _ (C6) 



Here, we consider the case with AQi > and Qi — Qj > as shown in Fig. IC1I The first two terms in equation (|C5[) . 
min{2(Qi — Qj), AQi}, ensure the condition of Qr > Qj (see Fig. IC1|I . This is the lower bound of Qr for AQi > 0. In 
the finite- volume method, the upper bound of Qr is determined by Q for s > Sj. On the other hand, in the SPH method, 
we do not know the distribution of Q for s > Si explicitly at the instance in the calculation of the interaction between i- 
and j-th particles. However, it can be estimated by the fact that dQi/ds is calculated by using all particles for s < Si and 
s > Si, suggesting that dQi/ds can be regarded as the average gradient around Xi. If the gradient in Sj < s < Si is simply 
approximated by (Qi — Qj)/Asij, the gradient in s > s it AQ'/Asy , can be guessed by equation (|C6|I . In this paper, we set 
the upper bound of Qr by AQ£ (see equation (IC5|) . 

In the actual calculation, we take into account the domain of dependence as follows, 

„ „ AQ, / C 3 At\ AQ, ( C l At\ 

In the RP, we use the monotonicity constraint with respect to p, v ■ n, P t , B\. In the MOC, the monotonicity constraint is 
used with respect to the Riemann invariants AJ± t , (see equation (|B4[) ). Using them, we derive AB±i and Av±.i. Equation 
(IC5|) suppresses numerical oscillations reasonably well. However, many improvements could still be made to the monotonicity 
constraint in the SPH method. 



